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Solutions of the semilinear wave equation are found numerically in three spatial dimensions with 
no assumed symmetry using distributed adaptive mesh refinement. The threshold of singularity 
formation is studied for the two cases in which the exponent of the nonlinear term is either p = 5 
or p — 7. Near the threshold of singularity formation, numerical solutions suggest an approach to 
self-similarity for the p — 7 case and an approach to a scale evolving static solution for p — 5. 
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I. INTRODUCTION 

One area of interest within the context of a nonlin- 
ear wave equation is the emergence of a singularity from 
smooth initial data. Among the issues raised by the for- 
mation of a singularity, the nature of the threshold for 
such formation is of particular interest. A number of 
past studies have addressed this threshold numerically 
in the nonlinear sigma model in both two 0, 0] and 
three 0, 0, IE IE spatial dimensions. Here, I study the 
semilinear wave equation following the work of || which 
considers the formation of singularities in the model re- 
stricted to spherical symmetry. 

A scalar field, <j), obeys the wave equation 
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(1) 



for p odd (preserving the symmetry <f> — > —(/>). In three 
dimensions with Cartesian coordinates, this equation be- 
comes 

4> = 4>,xx + 4>,yy + 4>,zz + <f> P , (2) 

where commas indicate partial derivatives with respect 
to subscripted coordinates and an ovcrdot indicates par- 
tial differentiation with respect to time. Solutions are 
found numerically by rewriting the equation of motion @ 
in first differential order form and replacing derivatives 
with second order accurate finite difference approxima- 
tions. These finite difference equations are solved with 
an iterative Crank-Nicholson scheme. In order to achieve 
the dynamic range and resolution needed to resolve the 
features occurring on such small scales in this model, 
adaptive mesh refinement is used. In fact, the code nec- 
essary for this is achieved with minimal modification to 
the code used in [E , and the reader is referred to these 
earlier papers for computational details. These changes 
consist of: (1) making the association <j> = x, (2) replac- 
ing the nonlinear term in the evolution update of that 
model with the last term of Eq. and (3) removing 
the regularity condition on the scalar field at the origin. 
As is common, an outgoing Robin boundary condition 
is applied at the outer boundaries assuming a spherical 
outgoing front. While not perfect, the outer boundary is 
generally far enough away from the central dynamics so 
that the boundary condition has no effect. 

Tests of the code include examining the convergence 
of various properties to the continuum properties as the 





Description 


4>(x,y,z,0) 


<j>(x,y,z,0) 


a 


Ellipsoid 


G 


or 

(yG, x -xG, y ) 


b 


Two pulses 


Gi + G-2 


SGi . dG 2 

Vl "a + V2~n 

ri.T * rltr. 


c 


Toroid 




{y<t>,m — x<t>,y) 


d 


Antisymmetric 


X G 





e 


Flat Pulse 


A (|tanh^ + i)x 
(itanh^ + I) 





f 


Static 
(for p = 5) 


G+(l+r 2 /3)- 1/2 






TABLE I: List of various initial data families. For fami- 
lies (a)-(f) both the field (j)(x, y, z, 0) and its time derivative 
4>(x, y, z, 0) are shown at the initial time t — in terms of var- 
ious parameters. The terms G, Gi, and G2 represent unique 
Gaussian pulses as defined in Eq. ©. In family (b), the pa- 
rameters vi and V2 are the respective velocities of the two 
pulses, generally chosen to have a grazing collision. Family 
(f ) represents perturbations of the static solution for the par- 
ticular case of p — 5. 



resolution is increased. To that end, the energy density 
associated with the scalar field is given by 



P 



4>) +{<t>, X f + (4>,yf + {4>,zf 



(3) 



so that the energy contained in the grid can be computed 
by an integral. Similarly, the z-componcnt of the angular 
momentum is H 



Jz 



d 3 x M xy = d 3 x <p {ycj) x - xepy) 



(4) 



One example of the conservation of these quantities is 
presented in Fig. Q Also in this figure is plotted the 
convergence factor 



Q(t) 



\Xih - X2h\2 
\X2h - Xh\2 



(5) 



which compares the differences in solutions as the reso- 
lution is increased. For a second order accurate scheme 
such as this one, the convergence factor is expected to 
converge to the value four. As shown in Fig. ^ the code 
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FIG. 1: Demonstration of convergence of the code for p = 7. 
The bottom frame displays the energy loss (relative to 
the initial time) AE(t) = E(t) — E(0) for four resolutions: 
(32 + l) 3 (solid, black crosses), (64+ l) 3 (dot, blue circles), 
(128 + l) 3 (short dash, green triangles), and (256 + l) 3 (long 
dash, red squares). As the resolution increases, energy loss 
decreases. In a similar fashion the middle frame shows the 
loss of the z-component of the angular momentum as a func- 
tion of time which also converges to conservation. The top 
frame shows the convergence factor Q(t) as in Eq.|S| which 
approaches the expected value of 4 with increasing resolution. 
This run is a bit below the threshold for singularity formation 
and, as such, represents a strong field example away from the 
linear regime. The initial data is from family (c) (see Ta- 
ble HJ with parameters A = 0.14, S = 2.5, e x = 2.4, e y = 0.4, 
f2 z = 3, and R = 3.8. The dimensionless ratio of the angular 
momentum to the energy squared is J/E 2 — 0.013. 

demonstrates second order convergence as well as conser- 
vation of the total energy and angular momentum. 

A number of families of initial data have been explored, 
and these are described in Table Q] Initial data is created 
by specifying <p and <j) & t the- initial time, and is done 
so here with a variety of real constants. A generalized 
Gaussian pulse is denned as 

G(x,y,z) = Ae~ {f - R)2/S \ (6) 
where f is a generalized radial coordinate 

r = \J £.t (x - x c f + e y (y - y c f + (z - z c ) 2 . (7) 

Such a pulse depends on parameters: amplitude A, shell 
radius i?, pulse width S, pulse center (x c , y c , z c ), and 



FIG. 2: Demonstration of threshold behavior of p — 5 static 
solution. The maximum value of the field is plotted for three 
evolutions; In solid line the results for the static solution 
(Family (f ) in Table |IJ are shown. In dashed line are the 
maximums obtained with a perturbation with negative am- 
plitude A = -0.05, R = 6, 5 = 2, e x = 2, and e y = 0.6. 
In dotted line the results with the same perturbation except 
with a positive amplitude A = +0.05 are shown. As the figure 
indicates, the perturbations send the solution towards either 
singularity formation or dispersal. Here, one can also see 
that at late times the unstable static solution is also driven 
to singularity formation because of the variety of numerical 
perturbations inherent in the numerical scheme such as the 
boundary treatment. 

skewing factors e x and e y . For e x ^ 1 +^ e y such a 
pulse has elliptic cross section. Family (a) represents 
a single pulse for which the parameter v takes the val- 
ues {— 1,0, +1} for an approximately out-going, time- 
symmetric, or approximately in-going pulse, respectively. 
The angular momentum of the pulse about the z-axis is 
proportional to the parameter Vt z as well as to (e x — e y ) 2 . 
The other families are similarly defined. 

II. THRESHOLD BEHAVIOR 

Solutions of Eq. approach one of two end states, 
either dispersal or singularity formation. This situation 
is quite similar to the evolution of a scalar wave pulse 
coupled to gravity which tends toward either dispersal 
or black hole formation. Choptuik's study of this prob- 
lem found fascinating behavior at the threshold for 
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FIG. 3: Example of a slightly super-critical evolution of an- 
tisymmetric initial data (Family (d) from Table for p — 5. 
The first frame shows the initial field configuration 4>{x, y, 0, 0) 
in the x-y plane and the second frame shows the two, well- 
resolved regions in which the field is blowing up. 



black hole formation (for reviews of black hole critical 
phenomena see [Til frj]'). It was largely in the spirit of 
Choptuik's work that the previously mentioned studies 
considered what happens at the threshold of singularity 
formation in the nonlinear sigma model PIS 0,00 SGI- 
Self-similar solutions are often found at the threshold, 
and indeed such is the case here. As mentioned in ||, 
the scaling symmetry of Eq. allows for self-similar 
solutions of the form 



cj ) (r,t) = (T-tr a U(p) 
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FIG. 4: Demonstration of the near critical approach to the 
static solution for p = 5. Shown are certain x > 0, y = 0, z — 
slices of an evolution at five different times near the collapse 
time, T. Each slice <f>(\rix) is rescaled to L 1 ' 2 c/>(\nx — InL) in 
accordance with the rescaling symmetry of the problem. The 
rescaling factor L(t) is chosen for each time slice in order to 
achieve unity at the origin. The (unrescaled) static solution 
us(lnr) is shown in solid line. The excellent agreement among 
these profiles suggests the solution represents a scale evolving 
static solution. The inset shows three different slices along the 
three positive axes all at the same time, and indicates that 
the central area of the solution is spherically symmetric. The 
initial data come from family (a) of Table[I]with R = 3, 8 = 2, 
x c = = y c = z c , e x = 2, e y = 0.7, v = 0, and Q z = 0.1, 
and the evolution entails 11 levels of 2:1 refinement with a 
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where 



T-t 



and 



p-1 



(9) 



Here, T is the collapse time associated with the forma- 
tion of a singularity. In 0, they find discrete families 
of solutions, U n (p), for p = 3 and for p = 7, but find 
no nontrivial self-similar solutions for p = 5. Plugging 
Eq. |(HJ| into Eq. 0), one arrives at an ODE that can 
be solved with a standard shooting method and which 
duplicates the results of @ . 

I now consider the different cases for p in turn. 



A. The case: p — 3 



that looked to be dispersing would eventually demon- 
strate growth near the origin, and hence the threshold of 
singularity formation could not be studied. Similar prob- 
lems are encountered here even without the assumption 
of spherical symmetry 



B. The case: p — 5 

In spherical symmetry 0] , no nontrivial self-similar so- 
lution exists for p = 5, and evolutions near the threshold 
approach the known static solution 



In the spherically symmetric evolutions of [8| for p = 3, 
the threshold could not be studied because evolutions 



(10) 
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FIG. 5: Near-critical solution from the explicitly spheri- 
cally symmetric code for p = 5 with initial data family 11 1 2t 
and JT3J) and R = 0.1, 5 = 0.3, and A* w 5.539. Similar to 
Fig. El the solution at times near the collapse time is shown 
rescaled. Its agreement with the (unrescaled) static solution 
(shown in solid line) suggest a scale evolving static solution. 



Retaining the notation of |8(, we have a family of static 
solutions, Ug, generated by rescalings of 1)10(1 



;(r) = L" 1 /V (r/L) 



(11) 



In three dimensions without the assumption of spheri- 
cal symmetry, the static solution remains at the threshold 
of singularity formation. In Fig. the maximum value 
of the scalar field within the computed domain is plotted 
versus time. For the static solution, this value remains 
essentially constant until late times when inherent nu- 
merical perturbations drive it to form a singularity. Also 
plotted are the results of initial data consisting of the 
static solution along with a small Gaussian pulse added 
explicitly to perturb the solution. For positive amplitude 
of the Gaussian pulse, the solution is driven to singular- 
ity formation whereas for negative amplitude the solu- 
tion disperses. These results suggest that the static solu- 
tion remains on threshold even without the assumption 
of spherical symmetry. 

However, for other families of initial data near the 
threshold, the evolution does not appear static. Instead, 
the collapsing region appears roughly self-similar in its 
collapse about some central point. Consider for example 
a family of initial data which is antisymmetric across the 
y-z plane, Family (d) from Table [IJ One such example is 



FIG. 6: Results from two near critical evolutions for p = 7 
with initial data family (e). The maximum of the field config- 
uration (f> is plotted versus time for a slightly subcritical evolu- 
tion (solid) and for a slightly super-critical evolution (dotted). 
The collapse time T is estimated by determining the lifetime 
of the supercritical evolution closest to criticality, in this case 
T m 3.67. The inset shows the same information with the 
expected scaling of <\>. 



shown in Fig. [3] The first frame in the figure shows the 
initial configuration for <j> and the second frame shows the 
solution near the collapse time. Two regions of collapse 
form with both regions spherically symmetric about their 
respective centers. 

That the collapse appears self similar when no such 
solution is admitted also occurs in the case of blowup 
with a Yang-Mills field 0, 0|. There, the dynamics 
were identified with a scale evolving static solution, and 
such an identification appears to be the case here. 

Evidence that the near critical solution represents the 
static solution with a scale factor dependent on time, 
L(t), is presented for a particular case in Fig. 0] Shown 
in the figure is the evolution at different times near the 
collapse time rescaled according to Eq. <|llfl by choosing 
L(t) such that the rescaled quantity VX0(r = 0) is unity. 
The (unrescaled) static solution, us, is also shown and 
the excellent agreement among the profiles is strong ev- 
idence that indeed the evolution is proceeding along the 
scale evolving static solution. 

The inset of Fig. 0] shows three different spatial slices 
of the data for a particular time. They agree quite well, 
providing good evidence that the solution is spherically 
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C. The case: p = 7 
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FIG. 7: Demonstration of self-similarity of a near-critical so- 
lution for p — 7. This is similar to Fig. [4] and indeed has 
the same initial data family (with different critical amplitude 
A*). That the solutions at different times t coincide indicates 
self-similarity. Also shown is the explicitly self-similar solu- 
tion Ui(p) found by solving an ODE. The collapse time T 
is determined by the average of the values required for each 
time slice to achieve agreement with Ui for p — 0. 



symmetric. 

One can observe similar behavior in spherical symme- 
try by modifying the code from Q . For many families of 
initial data, near threshold solutions approach the static 
solution in the conventional way. However, in Fig. [5] a 
near critical solution is shown obtained by tuning the 
initial data family 



#r,0) 
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tanh 
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1 



(12) 



with the initial time derivative set consistent with Eq. (jSJ) 



#r,0) , Mr,0) 

r T 

2 dr 



(13) 



For this family, near threshold solutions also appear to 
approach a scale evolving static solution. 



In [g] for the p = 7 case, they find in the critical limit 
an approach to the U\{p) self-similar solution. Here, 
without the assumption of spherical symmetry, similar 
threshold behavior is observed (see Fig.[SJ). Such a near 
critical solution is shown in Fig.[7| The solution appears 
both self-similar and spherically symmetric. The inset of 
Fig. compares the obtained solution to the ODE so- 
lution, Ui(p), and they appear quite similar suggesting 
that it remains the critical solution even without the as- 
sumption of spherical symmetry. 



III. CONCLUSIONS 

The semilinear wave equation represents perhaps the 
simplest nonlinear generalization of the linear wave equa- 
tion, and yet it displays interesting threshold behavior. 
In particular, this work has extended the results of Q 
obtained in spherical symmetry to the full 3D case. 

For the p = 3 case, the threshold could not be studied 
because of its late time growth which was also observed 
in spherical symmetry. 

For the p = 5 case, a critical solution is observed which 
appears roughly self similar despite the fact that no such 
solution is admitted. Instead, the solutions suggest an 
approach to the static solution with an evolving scale. 
This same behavior is also observed in a code which ex- 
plicitly assumes spherical symmetry. 

For the p = 7 case, a critical solution is found which 
resembles that found in the spherically symmetric case. 
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